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Abstract. We show that the widely used model governing the motion of two 
incompressible immiscible fluids in a possibly heterogeneous porous medium 
has a formal gradient flow structure. More precisely, the fluid composition is 
governed by the gradient flow of some non-smooth energy. Starting from this 
energy together with a dissipation potential, we recover the celebrated Darcy- 
Muskat law and the capillary pressure law governing the flow thanks to the 
principle of least action. Our interpretation does not require the introduction 
of any algebraic transformation like, e.g., the global pressure or the Kirchhoff 
transform, and can be transposed to the case of more phases. 


1. Introduction 

1.1. General motivations. The models for multiphase porous media flows have 
been widely studied in the last decades since they are of great interest in several 
fields of applications, like e.g. oil-engineering, carbon dioxide sequestration, or 
nuclear waste repository management. We refer to the monographs for an 

extensive discussion on the derivation of models for porous media flows, and to 
BIUDM for numerical and mathematical studies. 

More recently, F. Otto showed in his seminal work [T8] that the so-called porous 
medium equation: 

dtp — A p m = 0 for ( x , t) £ x R+ and m > 1, 

which is a very simplified model corresponding to the case of an isentropic gas 
flowing within a porous medium, can be reinterpreted in a physically relevant way 
as the gradient flow of the free energy with respect to some Wasserstein metric 
in the space of Borel probability measures. Extensions to more general degener¬ 
ate parabolic equations were then proposed for example in id ng. See also for 
instance [7j or jT5] for the interpretation of some dissipative systems as gradient 
flows in Wasserstein metrics. 

In this note, we will focus on the model governing the motion of an incompressible 
immiscible two-phase flow in a possibly heterogeneous porous medium, that will 
appear in the sequel as m and m da)- This model is relevant for instance for 
describing the flow of oil and water, whence the subscripts o and w appearing in 
the sequel of this note, within a rock that is possibly made of several rock-types. 
Our goal is to show that, at least formally, this model can be reinterpreted as the 
gradient flow of some singular energy. This will motivate the use of structure¬ 
preserving numerical methods inspired from [9] to this model in the future. 

Our approach is inspired from the one of A. Mielke m and, more closely, to 
the one of M. A. Peletier [T9]. The basic recipe for variational modeling is recalled 
in m then its ingredients are identified in £{2] This approach is purely formal, 
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but it can be made rigorous under some unphysical strict positivity assumption on 
the phase mobilities r] 0 ,r) w defined below. We will remain sloppy about regularity 
issues all along this note. 

1.2. The recipe of getting formal variational models. Here we recall very 
briefly the main ingredients needed for defining a formal gradient flow. 

i. The state space At is the set where the solution of the gradient flow can evolve. 

ii. At a point s £ At, the tangent space T s At, to whom would belong dts, is 
identified in a non-unique way with a so-called process space Z s (that might 
depend on s). More precisely, we assume that for each s £ At there exists an 
onto linear application V(s) : Z s — > T s At. 

iii. The energy functional £ : At — > R. U {+ 00 } admits a (local) sub-differential 
d s £(s) c (T S M)* at s £ At. 

iv. The dissipation potential T> is such that, for all s £ At and all ~V £ Z s , one 
has V(s; V) > 0. It is supposed to be convex and coercive w.r.t. to its second 
variable. 

v. The initial data s° belongs to At. 

All these ingredient being defined, we obtain from the principle of least action that 
s : R + -A At is the gradient flow of the energy £ for the dissipation V if 

(la) d t s = V(s)V 
where 

(lb) V € argmin ( max ( Visit)] V(t)) + (h . 7 ? (s)v\ 

\hed e £( s ) \ v ' \ /(t b m)*,t s m 

The formula m means that a gradient flow is lazy and smart: the motion aims 
to minimize the dissipation while maximizing the decay of the energy. We refer 
to m 119] for additional material on such a formal modeling and to [2J for an 
extensive (and rigorous) discussion on gradient flows in metric spaces. 



2. Variational modeling for two-phase flows in porous media 


2.1. State space and process space. Let f l be an open subset of representing 
a (possibly heterogeneous) porous medium , let <f> : fl — > (0,1) be a measurable 
function (called porosity ) such that (j) < (f>{x) < cfi for a.e. x £ fl for some constants 
(j>, (j) £ (0,1), and let s a ,s w : fl —> [0,1) be two measurable functions (so-called 
residual saturations) such that s a (x) + s w (x) < 1 for a.e. x £ fl. In what follows, 
we denote by 

s 0 (x) = 1 — s w (x), s w (x) = 1 — s Q (ai), for a.e. x £ £1. 


For almost all x £ fl, we denote by 


A™. = 


^ 'S — (s G , s w ) 


So + s w — 


1 with s a (x) < s a < s a (x) for a £ 



Let s° = (s[J,s°) be a given initial saturation profile, we denote by m a ( a £ 
{o, m}) the volume occupied by the phase a in the porous medium, i.e., 


m a = f (f(x)s 0 o (x)dx, and m w = f cfi(x)s 0 w {x)dx. 

Jn J q 
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For simplicity, we restrict our attention to the case where the volume of each phase 
is preserved: no source term and no-flux boundary conditions (otherwise, non- 
autonomous gradient flows should be considered). Hence the saturation profile lies 
at each time in the so-called state space A4, defined by 


M = { s = (s 0 ,s w ) 


: —K 


with / <f>(x)s a (x)dx = m a for a £ {o, w} 

Jn 


Let us now describe the processes that allow to transform the saturation profile. 
We denote by 



Vw) 


v a : fl —> 


with v a ■ n = 0 on 



the process space of the admissible processes for modifying a saturation profile 
s £ M. The identification between V = (v a , v m ) £ Z s and s = (s Q , s w ) £ T S M is 
made through the onto operator V(s) : Z s —>• T S M defined by 


( 2 ) 


V{s)\ 


-T V ' v o; -T V 



vv € z a . 


Since dts £ T S A4, the relation 0) yields the existence of some phase filtration 
speeds (v a ,v w ) £ Z s such that the following continuity equations hold: 


(3) 4>d t s a + V ■ v a = 0, a£{o,w}. 

The relation (0 must be understood as the local volume conservation of each phase 
a £ {o, w}. Finally, the duality bracket (•, ■) rr a M)*,T a M is given by 


(h, s) (T e M)* , TbM 




V h a • Voc 


2.2. About the energy. For a.e. x £ ft, we assume tv(-,x) : [s 0 (x), s 0 (*)] —> K 
to be a maximal monotone graph whose restriction 7T| (s _ (•, x) to the open interval 
(s 0 (x),s 0 (x)) is an increasing (single-valued) function belonging to L 1 (s 0 (x),s 0 (x)). 
In particular, 7r“ 1 (-,ai) : R —>• [s Q (ai), s D (ai)] is a single valued function. 

We denote by n : Rx —» RU{+oo} the (strictly convex w.r.t. its first variable) 

function defined by 

PS 0 

/ n(a,x)da - (p 0 - p w )sgz if s a £ \s 0 (x), s 0 (ai)], 

J cr(x) 

+oo otherwise, 


n(s D , x'j — 


where, denoting by e~ the downward unit normal vector of R w , we have set z = 
x e z , and where g and p a denote the gravity constant and the density of the phase 
a respectively, and where cr is such that x i —> ir(cr(x), x) — ( p a — p w )gz is constant. 
Since 7T| (a 3o) (-,aj) £ L 1 (s 0 (x),s 0 (x)), we get that n(s D (ai),ai) and n(s 0 (x),a;) are 
finite for a.e. x £ H. 

The volume energy function E : R 2 x fl —> R U {+oo} is defined by 


( 4 ) 


E(s, x) 


n(s G , x) if s (^o, Suf) £ , 

+oo otherwise. 
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The function E(-,x) is convex and finite on for a.e. x G £1. Its sub-differential 
is given by 


d s E(s,x) 


{(h 0 , h w ) g K 2 

0 


hw + (Po 


pw)gz G 7r 



if s G A x , 
otherwise. 


Finally, we can define the so-called global energy £ : Ad R U {+oo} by 


(5) £(s) = / (f(x)E(s(x), x)Ax, Vs = (s D , s w ) G Ad. 

J o 

The saturation profile s G .Ad is of finite energy £(s) < oo if and only if s( x) G 
for a.e. x G f2. For s G M with finite energy one can check that the local sub¬ 
differential d s £(s) off at s is given by 


(6) d s £(s) = | h = (h 0 , h w ) : fl —> R 2 s.t. 


h 0 — h w + (p 0 — p w )gz G 7r(s 0 , x) for a.e. x G f2 j. 


2.3. About the dissipation. The permeability tensor field A G L°°(£t-,M. NxN ) is 
assumed to be such that A(x) is a symmetric and positive matrix for a.e. x G $7. 
Moreover, we assume that there exist A*, A* G R+ such that 

A*|w| 2 < A (x)u ■ u < A*|it| 2 , for all u G R w and a.e. x G fi. 

This ensures that A(x) is invertible for a.e. x G ft. Its inverse is denoted by 

A _1 (a:). 

We also need the two Caratheodory functions g 0 ,Vw ■ R x —> R + — the 
so-called phase mobilities — such that g a (-,x) are Lipschitz continuous and non¬ 
decreasing on R + for a.e. x G fl and a G {o, w}. Moreover, we assume that 
7j a (s,x) = 0 if s < s a (x) and that g a (s,x) > 0 if s > s a (x). 

Given s = (s 0 ,s w ) G Ad and V = (v 0 ,v w ) G Z s , we define the dissipation 
potential T> by 

V(s,V) = ± J2 [ A 1V , a '” a ^ Vs G Ad, VV G Z s . 

2 a£rt Ja 

It is easy to check that dissipation is finite, i.e., T>(s, V) < oo, iff v a = 0 a.e. on {x G 
Q | s a (x ) < s Q (a:)}. 


2.4. Principle of least action and resulting equations. Let us consider the 
gradient flow governed by the energy £, the continuity equation Q, and the dissipa¬ 
tion T>. Let s G Ad be a finite energy saturation profile, then because of the principle 
of least action (00 and of the definition © of the operator V(s) : Z s — »• T s Ad, the 
process V = ( v a ,v w ) G Z s and the hydrostatic phase pressures h = ( h a ,h w ) must 
be chosen so that (V, h) is the min — max saddle-point of the functional 

(V,h)^V(s,V)- V [h a V-v a dx. 

Ar i An 


( 7 ) 
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One can first fix h G d s £(s) and minimize w.r.t. V. This provides 

( 8 ) 


Injecting this expression in 0 and maximizing w.r.t. h G d s £ (s), that is minimiz¬ 
ing 

(9) h = argmin f ^ f p a (s a )AVh a ■ Vh a 

had B £(s) V 

among all elements h in the subdifferential d s £(s), yields 

(10) - V • (v 0 + v w ^j = 0, v a =-p a (s a )AVh a . 

In (flOl) the first condition follows from the constraint h G d s £(s) in ©, and the 
second one from (JH]). 

Define the phase pressures p = (p 0 ,p w ) by p a {x) = h a (x) + p a gz , for a.e. x G D 
and a G {o, u;}, then we recover the classical Darcy-Muskat law: 

(11) v a = -?y Q (s Q )AV (p a - p a gz), a&{o,w}. 

Moreover, it follows from 0 that the following capillary pressure relation holds: 

(12) Po(x) — p w {x) G 7r(s 0 (x), x) a.e. in Q. 

We recover here the multivalued capillary pressure relation proposed in [EJ Uni- 
Combining (0 and (flOl) easily gives d t (s 0 + s w ) = 0, so that the condition 

(13) s 0 + s w = 1 a.e. in fi, 

is preserved along time and the whole pore volume remains saturated by the two 
fluids. 

Gathering ©, GH). CE2) and m gives the usual system of equations governing 
immiscible incompressible two-phase flows in porous media onus HUTU]. 




Remark 1. By similarity with the classical Wasserstein distance used in optimal 
mass transport 01 one could here endow the tangent space T s Ai at s G Ai with a 
weighted H^ 1 -scalar product 

(^i, S‘2^rp ^ ( I Pot (^a)-A- ‘ ^^2,ad3I, 

where, for i G {1,2} and a G {o,w}, we have set Si = (si,o, Si,™) and where hi^ a 
solves 


—V • (p a (s a )A'Vhi t0l ) = Si >a in fl, p a (s a )A'Vhi ta ■ n = 0 on dfl. 

Under some conditions on the functions p a (see [14]}, this should allow us to con¬ 
sider M as a metric space endowed with the corresponding distance, but £ is not 
locally X-convex for this Riemannian structure. The minimization © then consists 
in the selection of the subgradient with minimal norm. 
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